Pearson Chi-squared test with approximated methods for residual degrees of freedom in GLMMs

Explanations

package: glmmrBase https://cran.r-project.org/web/packages/glmmrBase/index.html

website (out of date): https://samuel-watson.github.io/glmmr-web/

preprint about the package (2024): https://arxiv.org/abs/2303.12657

mantainer: Samuel Watson - Prof stats Uni Birmingham https://www.birmingham.ac.uk/staff/profiles/applied-health/watson-sam

To create the Pearson Chi-squared test with the approximate residual degrees of freedom, I first:

  • converted the lme4 function to the glmmrBase function
  • and used the method small-sample_correction() to have the dfs correction:
    • Kenward-Roger
    • Kenward-Roger2
    • Satterthwaithe
  • I subtracted the sample size from this df to have the residual df.
  • I used two.sides test.

From Watson 2024: “The degrees of freedom correction is also given in Kenward and Roger (1997), which is the same degrees of freedom correction originally proposed by Satterthwaite (1946). The Kenward-Roger correction, though, can over-estimate the standard errors in very small samples in some cases. Kenward and Roger (2009) proposed an improved correction for covariance functions not linear in parameters (such as the autoregressive or exponential functions). The improved correction adds an additional adjustment factor. We have implemented these corrections in the package, which can returned directly from a Model object using the small_sample_correction() member function”

The function:

approximateDFpearson <- function(model, data, type=c("naive", "KR", "KR2", "sat"),
                              alternative=c("two.sided", "greater", "less")){
  
  if(class(model)[1] != "glmerMod") stop("model is not a glmerMod from lme4")
  if(model@resp$family$family != "poisson") stop("model is not a Poisson")
  
  if(type == "naive"){
    rdf <- df.residual(model) 
    } else{
    f1 <- glmmrBase::lme4_to_glmmr(formula(model), colnames(data))
    mod <- glmmrBase::Model$new(f1, data=data, family=poisson())
    rdf <- dim(data)[1]- mod$small_sample_correction(type)$dof[1]
    }
  
  rp <- residuals(model, "pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  
  if(alternative == "greater") pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  else if (alternative == "less") pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE)
  else if (alternative == "two.sided") pval <- min(min(pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE), pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)) * 2,1)
  
  out = list()
  out$statistic = prat
  names(out$statistic) = "dispersion"
  out$parameter = rdf
  names(out$parameter) = paste0("df-",type) 
  out$method = "Parametric dispersion test via mean Pearson-chisq statistic"
  out$alternative = alternative
  out$p.value = pval
  class(out) = "htest"
  return(out)
}
source(here("code_sim", "0_approximateDFpearson.R"))

Loading simulation outputs

Simulations done in 8_glmm_approxiDFpearson.R:

Sim setting:

sampleSize <- c(200, 500, 1000)
intercept <- c(-1.5,0,1.5)
overdispersion = seq(0,1,0.1)
ngroups <-  c(10, 50, 100)
nRep = 1000
slope = 1
randomEffectVariance = 1
load(here("data", "8_approximateDFpearson.Rdata")) # simulated data
(tspent <- map_dbl(out.pois, "time")) # minutes
##   10_200_-1.5      10_200_0    10_200_1.5   10_500_-1.5      10_500_0 
##      1.321605      1.219799      1.216456      2.353973      2.326321 
##    10_500_1.5  10_1000_-1.5     10_1000_0   10_1000_1.5   50_200_-1.5 
##      2.370980      8.441990      8.488686      8.568489      1.519092 
##      50_200_0    50_200_1.5   50_500_-1.5      50_500_0    50_500_1.5 
##      1.471159      1.405490      3.685771      6.909883      3.672749 
##  50_1000_-1.5     50_1000_0   50_1000_1.5  100_200_-1.5     100_200_0 
##     16.131561     16.184664     16.255666      1.805094      1.823181 
##   100_200_1.5  100_500_-1.5     100_500_0   100_500_1.5 100_1000_-1.5 
##      1.689025      5.617923      5.463803      5.451093     28.755045 
##    100_1000_0  100_1000_1.5 
##     31.243332     28.750171
#simulations
simuls.pois <- map_dfr(out.pois, "simulations", .id="model") %>%
  separate(model, c("ngroups", "sampleSize", "intercept"), sep="_") %>%
  rename(overdispersion = controlValues,
        PearNaive.rdf = PearNaive.rdf.df.naive,
        PearSat.rdf = PearSat.rdf.df.sat,
        PearKR.rdf = PearKR.rdf.df.KR,
        PearKR2.rdf = PearKR2.rdf.df.KR2)

There were many negative residual dfs, it generates NaN to the test:

simuls.nan <- simuls.pois %>% filter(PearSat.rdf < 0 |
                                     PearKR.rdf  < 0 |
                                     PearKR2.rdf < 0)

It was around 6% of the simulations.

Almost everything came from the 100 groups in 200 sample sizes - 2 observations per group:

table(simuls.nan$ngroups, simuls.nan$sampleSize)
##      
##         200   500
##   100 15567   392
##   50   1572     0

The percentage of failure was around the same among the corrections

sum(simuls.nan$PearSat.rdf <0) / dim(simuls.pois)[1]
## [1] 0.02415152
sum(simuls.nan$PearKR.rdf <0) / dim(simuls.pois)[1]
## [1] 0.02332997
sum(simuls.nan$PearKR2.rdf <0) / dim(simuls.pois)[1]
## [1] 0.02317508

Subsetting the simulation output only for the positive rdfs:

simuls.pois2 <- simuls.pois %>% filter(PearSat.rdf > 0 &
                                                       PearKR.rdf >0  &
                                                       PearKR2.rdf >0)

Comparing corrections of residual df

dof <- simuls.pois2 %>% dplyr::select(contains("rdf"), replicate, ngroups,
                                        overdispersion, intercept, sampleSize) %>%
  pivot_longer(1:4, names_to = "test", values_to = "res.df") %>%
  group_by(sampleSize, ngroups, intercept,overdispersion, test) 
dof$intercept <- fct_relevel(dof$intercept, "-1.5", "0", "1.5")
dof$ngroups <- fct_relevel(dof$ngroups, "10", "50", "100")
dof$sampleSize <- as.factor(as.numeric(dof$sampleSize))

Very similar rdfs. Satterthwaite had a bit smaller rdfs for small sample size.

Type I error

p.pois <- simuls.pois2 %>% dplyr::select(ends_with(".p"), replicate,
                                        ngroups,
                                        overdispersion, intercept, sampleSize) %>%
  pivot_longer(1:4, names_to = "test", values_to = "p.val") %>%
  group_by(sampleSize, ngroups, intercept, overdispersion, test) %>%
  summarise(p.sig = sum(p.val<0.05,na.rm=T),
            nsim = length(p.val[!is.na(p.val)]))
## `summarise()` has grouped output by 'sampleSize', 'ngroups', 'intercept',
## 'overdispersion'. You can override using the `.groups` argument.
p.pois$prop.sig <- p.pois$p.sig/p.pois$nsim
p.pois$intercept <- fct_relevel(p.pois$intercept, "-1.5", "0", "1.5")
p.pois$ngroups <- fct_relevel(p.pois$ngroups, "10", "50", "100")
p.pois$sampleSize <- as.factor(as.numeric(p.pois$sampleSize))

Power

  • Power for the 3 correction methods is basically the same.
  • Corrected power still very bad for small intercepts (-1.5 and 0),

figure for each N groups separated

Dispersion statistics

Median of disp stat

d.pois <- simuls.pois2 %>% dplyr::select(ends_with("dispersion"), replicate, ngroups,
                                        overdispersion, intercept, sampleSize) %>%
  pivot_longer(1:4, names_to = "test", values_to = "dispersion") %>%
  group_by(sampleSize, ngroups, intercept,overdispersion, test) %>%
  summarise(median.stat = median(dispersion, na.rm=T))
## `summarise()` has grouped output by 'sampleSize', 'ngroups', 'intercept',
## 'overdispersion'. You can override using the `.groups` argument.
d.pois$intercept <- fct_relevel(d.pois$intercept, "-1.5", "0", "1.5")
d.pois$ngroups <- fct_relevel(d.pois$ngroups, "10", "50", "100")
d.pois$sampleSize <- as.factor(as.numeric(d.pois$sampleSize))

Simplified figure with just one correction (KR2) to compare with the Naïve

Compare power with the other tests

Compare pearson corrections with Dharma-conditional and pearson par boot.

Using only ONE correction, 100 groups

load(here("data", "5_glmmPois_power_50.Rdata")) # simulated data
out.pois50 <- out.pois
load(here("data", "5_glmmPois_power_100.Rdata")) # simulated data
out.pois100 <- out.pois

out.pois <- flatten(list(out.pois50, out.pois100))

simuls.pois <- map_dfr(out.pois, "simulations", .id="model") %>%
  separate(model, c("ngroups", "sampleSize", "intercept"), sep="_") %>%
  rename("overdispersion" = "controlValues")

# power
pw <- simuls.pois %>% dplyr::select(dhaCO.p.val,refCO.p.val, replicate,
                                      ngroups,
                                      overdispersion, intercept, sampleSize) %>%
  filter(intercept %in% c( "-1.5", "0", "1.5"),
         sampleSize %in% c("200", "500", "1000")) %>%
  pivot_longer(1:2, names_to = "test", values_to = "p.val") %>%
  group_by(sampleSize, ngroups, intercept, overdispersion, test) %>%
  summarise(p.sig = sum(p.val<0.05,na.rm=T),
            nsim = length(p.val[!is.na(p.val)]))
## `summarise()` has grouped output by 'sampleSize', 'ngroups', 'intercept',
## 'overdispersion'. You can override using the `.groups` argument.
pw$prop.sig <- pw$p.sig/pw$nsim
pw$intercept <- fct_relevel(pw$intercept, "-1.5", "0", "1.5")
pw$sampleSize <- as.factor(as.numeric(pw$sampleSize))
comb <- p.pois %>% filter(ngroups %in% c("50", "100"), test == "PearKR2.p") %>%
  bind_rows(pw)
  
ggplot(comb, aes(x=overdispersion, y=prop.sig, col=test, linetype=ngroups, shape=ngroups))+
  geom_point(alpha=0.7) + geom_line(alpha=0.7) +
  facet_grid(sampleSize~intercept) +
  geom_hline(yintercept = c(0.05,0.5), linetype="dotted") +
  ylim(0,1)+
  theme(panel.background = element_rect(color="black"),
        legend.position = "bottom") + 
  guides(color=guide_legend(nrow=3, byrow=TRUE),
         linetype=guide_legend(nrow=2, byrow=TRUE),
         shape=guide_legend(nrow=2, byrow=TRUE)) +
  xlab("Overdispersion") + ylab("Power")+
  scale_color_discrete(labels = c("Sim-based conditional", 
                                   "Pearson Chi-sq KR2", 
                                   "Pearson Param. Boot."))

ggsave(here("figures", "8_glmm_approxDF_power2.jpeg"), width=8, height = 8)